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We determine optimal designs for some regression models which 
are frequently used for describing three-dimensional shapes. These 
models are based on a Fourier expansion of a function defined on 
the unit sphere in terms of spherical harmonic basis functions. In 
particular, it is demonstrated that the uniform distribution on the 
sphere is optimal with respect to all <3? p criteria proposed by Kiefer in 
1974 and also optimal with respect to a criterion which maximizes a p 
mean of the r smallest eigenvalues of the variance-covariance matrix. 
This criterion is related to principal component analysis, which is the 
common tool for analyzing this type of image data. Moreover, discrete 
designs on the sphere are derived, which yield the same information 
matrix in the spherical harmonic regression model as the uniform 
distribution and are therefore directly implementable in practice. It 
is demonstrated that the new designs are substantially more efficient 
than the commonly used designs in three-dimensional shape analysis. 

1. Introduction. Over the last decade, tools for acquiring and visual- 
izing three-dimensional (3D) models have become integral components of 
data processing in many fields, including medicine, chemistry, architecture, 
agriculture and biology. Volumetric shape analysis permits an evalutation 
of the actual structures that are implicitly represented in 3D image data. 
For the analysis, description and comparison of shapes of various struc- 
tures, shape descriptors, which are able to handle very different shapes and 
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to represent their global and local features, are of increasing interest (see, 
e.g., Brechbuhler, Gerig and Kiibler [4], Novotni and Klein [14], Szekely, 
Kelemen, Brechbuhler and Gerig [23], Ding, Nesumi, Takano and Ukai [5], 
Funkhouser, Min, Kazhdan, Chen, Halderman, Dobkin and Jacobs [7] and 
Kazhdan, Funkhouser and Rusinkiewicz [11], among many others). Spher- 
ical harmonic shape descriptors usually describe the surface in terms of a 
relatively small number of coefficients of a spherical harmonic expansion of 
the radius as a function on the unit sphere (see, e.g., [5] or [4]), that is, 

oo t 

(1.1) r(M) = E E « m (M), 

£=0 m=-l 

where 6 6 [0, 7r], £ (— vr,7r], the quantities 

(1.2) cf = — r T r(6,<p)Y e m (9,4>)d(t)sm6d6 

4tt Jo 

are the usual "Fourier" coefficients and 

{Y e m (6, <t>)\m e {-£, -£+l,...,£};£e N } 

is a complete orthonormal basis on the unit sphere. Let = r(6i, (pi) denote 
the observed radius of the 3D shape at polar angle 9{ and azimuthal an- 
gle (pi [in other words, the corresponding point of the shape has spherical 
coordinates (rjsin^jcos^rjsin^sin^jjrjcos^i) 71 ] and assume that data 

{{ri,9i,(j)i)\i = l,...,n} 

are available for one object. Usually a truncated expansion of order d is ap- 
plied as an approximation of (1.1), where the coefficients c™ are determined 
by the least squares critierion 

{n / d £ 

Eh-E E <f Y n0iAi) 
i=l \ e=0m=-£ 

and the estimated coefficients in this expansion (appropriately normalized) 
are then used for describing and analyzing the 3D shapes. For this purpose 
typical tools from multivariate statistics (cluster, discriminant or principal 
component analysis) are applied (see Ding, Nesumi, Takano and Ukai [5], 
Kazhdan, Funkhouser and Rusinkiewicz [11] and Kelemen, Szekely and 
Gerig [12], among many others). 

A common experimental design for this situation is a uniform distribution 
on the rectangle [0,7r] x (— tt,tt] realized by a grid or a uniform design that 
takes observations on several circles with equal distance on the z axis (see, 
e.g., [5]). In the literature on shape analysis these designs are mainly moti- 
vated by their easy implementation. If the grid is fine enough or a sufficiently 
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large number of circles on the unit sphere are used, the matrix B B of the 
least squares estimate c = {B T B)~ l B T r approximates a diagonal matrix, 
which simplifies the numerical calculation in the statistical analysis. Here 
r = (ri, . . . ,r n ) is the vector of measured radii and 

B = (Y e m (9 i Ai))f m 

is the design matrix corresponding to the least squares problem (1.3) (see 
[4])- 

In the present paper we consider the problem of finding optimal designs 
for 3D shape analysis based on spherical harmonic descriptors. In Section 2 
we present some more details on spherical harmonic descriptors and basic 
results on the theory of optimal experimental design. In the same section 
we also demonstrate that the uniform distribution on the unit sphere is 
optimal with respect to any of Kiefer's [13] $ p criteria if the interest of the 
experimenter is the estimation of the complete vector 

c=(c° ,c^A,cl...,c- d d fGR^ 2 

or a certain subset of the parameters. It is also shown that for any t < (d+l) 2 
this design maximizes the pth. mean of the r largest eigenvalues of the 
variance-covariance matrix {B T B) -1 . Therefore the uniform distribution on 
the sphere is particularly efficient for principal component analysis, which 
is the main tool for summarizing the information contained in the spheri- 
cal harmonic coefficients obtained by (1.3) (see, e.g., [5] or [11]). Because 
this design is continuous with density j^smO dOdcj), it is not directly im- 
plementable in practice. Therefore, for a finite sample size we determine in 
Section 3 discrete designs which give the same information matrix as the 
uniform distribution on the sphere. For this reason these designs are also 
optimal with respect to Kiefer's [13] & p criteria, and optimal with respect 
to a criterion that maximizes a p mean of the r smallest eigenvalues of the 
information matrix and is related to principal component analysis. In Sec- 
tion 4 we present several examples which illustrate the advantages of our 
approach, and determine optimal uniform designs which take in each direc- 
tion only one observation and are for this reason particularly attractive to 
practitioners. We also reanalyze a design used by Ding, Nesumi, Takano and 
Ukai [5] for principal component analysis and demonstrate that the designs 
derived in the present paper are substantially more efficient. Finally, some 
conclusions and directions for future research are mentioned in Section 5, 
while more technical details are given in the Appendix. 

2. Spherical harmonic descriptors and optimal design. An orthogonal 
system 



{Y e m (d, (/>)]£ G No; m e {-£, -£+!,...,£}} 
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of functions on the unit sphere satisfies 



(2.1) -!- r r Yp{e,(j))Yf '(0,<j>) #sin0 dO = 

4lT Jo J -it 

whenever (m, 7^ (m', £') . The common system used in shape analysis (see [4] 
or [5]) is obtained by the normalization 



* (Yf 1 (6, (j>)fd<t> siu6 dO = l 



1 

47T Jo 

and given by the spherical harmonic functions 

Y°(6,4>) = V / 2nH r lP°(cos0), n G N 



Y?(e,<f>) = </2(2n + 1)^— ^P™(cos0)cos(m<£), 
y {n + m)\ 

(2.2) m = 0,...,n;neN, 



Y" m (M) = J2(2n + l) i n + m || Pn m (cosg)sin(m0), 
y (n — mj! 

m = —n, . . . , — 1; n E N, 

where P™{x) is the mth associated Legendre function of degree n satisfying 
the differential equation 

(2.3) (1 - x 2 )P"(x) - 2xP'{x) + |n(n + 1) - ^ ™ ^ |p(g) = 0. 

It is well known (see [1], Chapter 9) that these functions can be repre- 
sented as 

(2.4) p-(x) = (-ir|^(i - x'r /2 ct + m 1/2 \x), 

2rm\ 

where cj?\x) is the kth ultraspherical polynomial orthogonal with respect 
to the measure (1 — x 2 ) a ~ 1 / 2 dx (see [22]). Note also that P°(a;) is the nth 
Legendre polynomial P n (x) orthogonal with respect to the Lebesgue measure 
on the interval [—1,1]. Because orthogonal polynomials can be calculated 
recursively, this representation allows a fast computation of the functions 
Y™, and the first four spherical harmonic functions corresponding to the 
case d = 1 are given by 

y o °(0,0) = i, Yi(e,<t>) = V3c OS e, 

(2.5) 

Y" 1 (#,(/>) = V3 sin 6 sin^, Y^{Q,<f>) = V3sm9cos(j>. 
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Figures 1 and 2 show the spherical harmonic descriptors 

(Y e m {9, <f>) sin 9 cos 0, Y[ n {9, <f>) sin 9 sin </>, Y e m (9, <p) cos 9) T ', m = -£,...,£, 

for ^ = 0, 1, 2, 3, when (0, </>) varies in the rectangle [0, ir] x (— 7r, tt]. 

Consider the regression model corresponding to the least squares prob- 
lem (1.3), 

(2.6) E[Y\9^]=c T f d {9^), Var[Y|M = a 2 > 0, 



G 
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where 

/^^) = (y o °(e^),yf 1 (0,0),y 1 °(0,0), 

(2.7) 

Yl (e, 0) , . . . , Y d ~ d (e, 0), . . . , Y d d (e, e k^ 1 ) 2 

is the vector of spherical harmonic functions of order d and 
c=(clci\c° 1 ,c\,...,c- d d ,...,c d d feR^ 2 

is the corresponding vector of parameters. Note that (d + l) 2 spherical har- 
monic functions appear in the model (2.6). An approximate design is a prob- 
ability measure on the set [0,tt] x (— w,tt]. For a probability measure with 
finite support, the support points, say Z{ = (8i,cf>i), determine the points on 
the sphere where the radius of the 3D shape is observed, and the correspond- 
ing weights, say Wi, give the relative proportion of total observations taken 
in a particular direction. For a given design £ with finite support, the covari- 
ance matrix of the least squares estimate for the vector c is approximately 
proportional to the information matrix 

(2.8) M(0 = r r f(9, <P)f T (9, 0) #(0, 4>) 

J-n JO 
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(this is essentially the matrix B B mentioned in the Introduction), and 
an optimal approximate design maximizes an appropriate function of this 
matrix. There are numerous criteria proposed in the literature that can be 
used to discriminate between competing designs (see [19] or [16]), and we 
will restrict ourselves to the famous family of & p criteria introduced by 
Kiefer [13], and a new optimality criterion which has not been considered so 
far in the literature and is related to principal component analysis. 

To be precise, let KeM.^ d+l ^ xs denote a given matrix of rank s < (d + 1) 
and assume that the main interest of the experimenter is in the estimation of 
s linear combinations K T c. If n observations are taken according to an ap- 
proximate design (possibly by applying an appropriate rounding procedure 
to a design with finite support or by using a discretization of a continuous 
design; see [17]), then the covariance matrix of the least squares estimate 
for K T c is approximately given by 

2 

-{K T M~{OK), 
n 

where A~ denotes a generalized inverse of the matrix A and we assume that 
the linear combinations K T c are estimable by the design £, that is, 

range(iT) C range(M(£)); 

see [16]. Let — oo <p < 1. Following [13], we call the design £* ^-optimal 
for estimating the linear combinations K T c if £* maximizes the expression 

(2.9) %(C) = (tr(K T Ar(OK)- p f p 

among all designs for which K T c is estimable. If K = /( rf+1 )2 is the identity 
matrix of order (d + l) 2 x (d + l) 2 , then £* is briefly called <J> p -optimal. Note 
that the cases p = and p = — oo correspond to the frequently used D- and 
S-optimality criteria, that is, 

$ (0 = det(K T M-(OK)-\ = X^d^M^iOK)- 1 ), 

(2.10) 

while the A criterion is obtained for the choice p = — 1, that is, 
(2.11) $_ 1 (£)=tr(K T M-(Z)K)- 1 . 

In the following discussion, we are interested in a design that is particularly 
efficient for the estimation of the coefficients corresponding to the (2k + 1) 
spherical harmonic functions 

V-k v -k+l v k 
Y k 

of the vector of regression functions defined in (2.7), where k G {0, ...,d} 
denotes a given "level of resolution." For this, define 0fc iS as the (2k + 1) x 
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(2s + 1) matrix with all entries equal to 0, let q G No, consider the indices 
< ko < k\ < &2 < • • • < k q < d and define the matrix 

(2-i2) h 1 hi iA7,)j 

by 

(2.13) Ki^i **' V-^ k u 

(i = 0, . . . , q; j = 0, . . . , d) . Note that ifGR (d+1)2xs , where s = £ | =0 (2A* + 1) , 
and that K T c gives the vector with coefficients 

(2.14) {cf e \m e {-k e , -k e + l,..., kt};£= 0, . . . ,q}. 

Two cases are of particular interest here and are therefore mentioned sepa- 
rately. If q = d, then 

(2-15) Kl^ d = I (d+1)2 

and precise estimation of the full vector of parameters is the main goal for 
the construction of the optimal design. On the other hand, if q = 0, only 
the coefficients corresponding to the (2ko + 1) spherical harmonic functions 
Y^~ fc ° , . . . , are of interest and the corresponding matrix is given by 

(2.16) /^ = [0 fe0i0 : ••• \Q koM . 1 \h M : ••• lO^jei^ 1 ^^ 2 . 

Note that the general matrix defined by (2.12) consists of (g + 1) (block) rows 
of the form (2.16). The following result shows that for matrices of this form, 
the <l?p criteria defined in (2.9) are maximized by the uniform distribution 
on the sphere independently of p £ [— oo, 1). 

Theorem 2.1. Let p G [— oo, 1), let < ko < ■ ■ ■ < k q < d be given indices 
and denote K = iffo,...^ as the matrix defined by (2.12) and (2.13). A <5 p - 
optimal design for estimating the linear combinations K T c in the spherical 
harmonic regression model (2.6) is given by the uniform distribution on the 
sphere, that is, 

(2.17) C{d0,d(t}) = — sm 9d9d(f). 

Moreover, the corresponding information matrix in the spherical harmonic 
regression model is given by M(£*) =^(2d+i) 2 - 

Proof. Let £* denote the design corresponding to the density (2.17). 
Then, due to the orthonormality of the spherical harmonic functions Y™, it 
follows that 

M(f ) = {[ f Y e m (e,4>)Yp , (e,4>)d ( p S m9^-) =/ {M+1)2 . 



1,1' ,m' ,m 
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Assume for the moment that p > — oo. According to Theorem 7.20 in [16], 
we obtain that the measure £* is <3? p -optimal if and only if the inequality 

(2.18) fJ(e^)K(K T K)-P- l K T f(e,^<tx(K T K)-P 

holds for all 9 G [0,7r] and <fi G (— tt, 7r]. Now the special structure of the 
matrix K defined in (2.12) and (2.13) yields 

K T K = I S , 

where s = X^=o(2^ + !)• Observing the definition of the vector fd(9,4>) 
in (2.7) and again the definition of the matrix K, (2.18) reduces to 

*>E E (W<«) 2 

£=0 m=—ki 

(2-19) 

= ±(2k e + 1){(P,W)) 2 + 2 £ 7TT^(P^os9)A 

for all 9 G [0, tt] and <f> G (— 7r, 7r], where we have used the representation (2.2) 
and the trigonometric identity cos 2 (mc/)) + sin 2 (m4>) = 1 for the last equality. 
From the identity 

P® (cos a cos P + sin a sin /3 cos <^>) 

= P,°(cosa)P fc °(cos/3) + 2 E ^^Pr(™sa)Pf (cos /3) cos(m^ 

m=l ^ ' 

for the Legendre functions (see [1], page 457) and the fact that P® is the 
kth. Legendre polynomial, it follows (with the choice a = [3, cj) = 0) that 

(P fe °(cosa)) 2 + 2 E f^(^(cosa)) 2 = P fe °(l) = 1. 

771=1 >• ' 

Now the right-hand side of (2.19) simplifies to 

E(2fc£ + l) = s 

£=0 

and, consequently, (2.19) or equivalently (2.18) holds for all 9 G [0,7r] and 
cf) G (— 7r,7r]. This proves that the design corresponding to the uniform 
distribution on the sphere is <E> p -optimal for any p G (— oo, 1) and any ma- 
trix K of the form (2.12) and (2.13). The remaining case p = — oo finally 
follows from Lemma 8.15 in [16]. □ 

It should be noted at this point that the E criterion is of particular im- 
portance for the design of experiments in 3D shape analysis, because many 
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authors propose to use principal component analysis for the comparison of 
different shapes (see, e.g., Ding, Nesumi, Takano and Ukai [5] or Kazhdan, 
Funkhouser and Rusinkiewicz [11], among others). More precisely, it is pro- 
posed to summarize the information contained in the spherical harmonic co- 
efficients by a principal component analysis based on the variance-covariance 
matrix of the least squares estimator obtained from (1.3), which explains 
the 3D shape variation of different objects by the first, say r, principal com- 
ponents corresponding to the r largest eigenvalues of the matrix M -1 (£). 
Consequently, an efficient design for principal component analysis should 
minimize a function of the largest r eigenvalues of the matrix M" 1 ^). To 
be precise, let r < (d + l) 2 and — p < oo < 1. Then we call a design £* ^f p<r - 
optimal if £* maximizes 

/ r \ VP 

(2.20) *pAo=[E{ x u)( M m p ) » 

where A(j)(M(£)) denotes the jth smallest eigenvalue of the matrix M(£*). 
Again the cases p = — oo and p = are obtained from the corresponding 
limits, that is, *_ oo , r (e) = A ( i ) (M(0) and tf 0|t .(£) = U r j= i A (j) (M(£)), re- 
spectively. To our knowledge, optimality criteria of this type have not been 
considered in the literature so far. Note also that the E- and A-optimality 
criteria (with K = I^+i) 2 ) are obtained for the choices r = 1 and r = (d+ l) 2 
with p = — 1, respectively. The most important case for the choice of the 
parameter p£ [— 00, 1) is certainly obtained for p= —1. In the following 
corollary, we show that the uniform distribution on the sphere is also ^ r ,p- 
optimal. 

Corollary 2.2. Let p e [—00, 1) and 1 < r < (d + l) 2 . Then the uni- 
form distribution on the sphere defined by (2.17) is ^ p ^-optimal in the spher- 
ical harmonic regression model (2.6). 

Proof. The case p = — 00 is obtained from Theorem 2.1 with K = 
I( d+ i)2. For -00 <p< 1, note that % >r (C) = rl/p if p/0 and *o,r(£*) = 1- 
If p ^ and were not $ p r -optrmal, then there would exist a design, say 
£, with 

(%A0Y = E A£,(M (fl) < r = (%, r (e)) p - 

3=1 

Consequently, we obtain Am(M(£)) < 1 = A m ; n (M(^*)), which implies that 
the uniform distribution on the sphere £* is not S-optimal and contradicts 
Theorem 2.1. The case p = is obtained by a similar argument and the 
assertion of the corollary has been established. □ 
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3. Discrete optimal designs. Note that the uniform distribution defined 
by (2.17) is not directly implementable as a design in real experiments. 
Therefore, for practical applications it is important to obtain discrete 
designs £ which are equivalent to the uniform distribution £* (d9 , d(f>) = 

sin dO dcj) in the sense that 

(3-1) M(0 = M(r) = / (d+ i) 2 , 

where M(£) is the information matrix (2.8) of the design £ in the spherical 
harmonic regression model (2.6). Note that due to Caratheodory's theorem 
(see [19]), there always exist discrete designs satisfying (3.1), and in the 
following we will construct a broad class of discrete designs that can easily 
be implemented in practice. For this purpose we need an auxiliary result 
about quadrature formulas, which will be proved in the Appendix. 

Lemma 3.1. Assume that r G N and let —1 < x\ < X2 < ■ ■ ■ < x r < 1 
denote r points with corresponding weights wi,...,w r > (J2i=i w i = !)• 
The points Xi and weights Wi generate a quadrature formula of degree z>r, 
that is, 

r r 1 

(3.2) J2wixf = ^ x e dx, £ = 0,...,z, 

if and only if the following two conditions are satisfied: 

(A) The polynomial V r (x) = ni=o( 2 ' — x i) ^ s orthogonal to all polynomials 
of degree z — r with respect to the Lebesgue measure, that is, 

(3.3) / V r (x)x e dx = 0, £ = 0,...,z-r. 



-i 

(B) The weights Wj are given by 
(3.4) Wj = \ f £j(x)dx, 



where 



= n 



denotes the ith Lagrange interpolation polynomial with nodes Xi,...,x r . 

It follows from condition (3.3) that z < 2r — 1 (otherwise there does not 
exist a solution of this system). Moreover, there exists at least one set of 
points that satisfies condition (3.3), namely 

(3.5) {x\P r (x)=0}, 
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where P r denotes the rth Legrendre polynomial orthogonal with respect to 
the Lebesgue measure on the interval [—1,1], because it is well known that 
this polynomial has r distinct roots located in the interval (—1, 1) (see [22]). 
Moreover, the corresponding weights defined by (3.4) are positive and define 
a Gaussian quadrature formula, that is, a quadrature formula of degree 2r — 1 
with r nodes (see, e.g., [6] or [8]). In the following discussion we will use 
this lemma with z = 2d. Then it follows that for any r G {d + 1, . . . , 2d}, 
there exists at least one quadrature formula {xi,Wi}l =1 determined by the 
equations (3.3) and (3.4) that integrates polynomials of degree 2d exactly 
[in other words, the system of equations in (3.2) is satisfied with z = 2d] 
and the corresponding weights Wi are positive. We consider a quadrature 
formula of this type. Define 



(3.6) #i = arccos£i, i = l 
and consider the design 

(3.7) fj, 



ui ■ ■ ■ u r 
W\ ■ ■ ■ w r 



on the interval [0, it]. Similarly, we define for any t E N and any a E (— ^=7r, — n] 
a design v = u(a,t) on the interval (— 7r,7r] by 




(3.8) 

where the points <j)j are given by 



t 



2tti 

(3.9) </,. = (*+—£, j = l,...,t. 

The following result shows that designs of the form \i (g> v are discrete & p - 
and ^> p ^-optimal designs for the spherical harmonic regression model (2.6). 



Theorem 3.2. Let p e [-co, 1), let < k < ■ ■ ■ < k q < d and denote by 
K = Kkg^k the matrix defined by (2.12). For any t>2d+l and any r G 
{d+ 1, . . . , 2d}, the design n®v with factors given by (3.7) (corresponding to 
a quadrature formula of degree 2d) and (3.8) is Q p - optimal for estimating the 
coefficients K T c, and is ^ p ^ r -optimal in the spherical harmonic regression 
model (2.6). 

Proof. Observing the proofs of Theorem 2.1 and Corollary 2.2, the 
assertion can be established by showing the identity 



(3.10) 



M(/i<g>z/) = / (d+1)2 . 
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For this let 

rp 

i>(<t>) = {^~d{(t>)^-d+\{(t>),---^d{4>)) 

= (V2 sm(d(j)) , . . . , V2smc/), 1, \/2cos</>, . . . , \/2 cos(d(f>)) T . 

Then the regression functions in the spherical harmonic regression model of 
degree d are given by 

(3.11) 7i ^(cos0)-^-(0), j = -i,...,i,i = 0,...,d, 

where P- is the Legendre function defined by (2.4) and the constants 7^ are 
given by 



7ii = 4/(2f + l) (l 



[Note that the different scaling of the cases j = and j 7^ in (2.2) has 
been accommodated by introducing the factor v2 in the definition of the 
functions iftj.] Therefore, the identity (3.10) is equivalent to the system of 
equations 

(3.12) lijlM f r Pl jl (cos 6)^ (<j))Pl el (cos e)Tp t (<j>)dti(0)du( ( j>) = 8 ik 6 je 

J~w JO 

(j = —i, . . . = 0, . . . , d;£ = —k, . . . ,k;k = 0, . . . , d), where 5ik denotes 
Kronecker's symbol. Observing Fubini's theorem, this system is satisfied 
if the equations 

(3.13) 7ijlM P? (cos0)P|(cos0) dfi(6) = 5 ik 5 je 

Jo 

(j = 0, . . . , i; i = 0, . . . , d; I = 0, . . . , k; k = 0, . . . , d) and 

(3.14) r^j(MM0)<M$=*i* j,£ = -d,...,d, 



can be established. It is well known (see Pukelsheim [16]) that the last 
identity is satisfied for measures of the form (3.8). To prove the remain- 
ing identity (3.13), recall that the measure [i corresponds to a quadrature 
formula that integrates polynomials of degree 2d exactly. Moreover, it fol- 
lows from the representation (2.4) that for an even m E {0,...,n}, the 
Legendre function P™(x) is a polynomial of degree n, while for any odd 
m G {0, . . . , n}, the function P™(x) / \J\ — x 2 is a polynomial of degree n — 1. 
Observing the equation (3.14), we can restrict ourselves to the case j = i for 
which the integrand P- (x)P^(x) in (3.13) is always a polynomial of degree 
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i + k < 2d, which contains the factor (1 — x 2 ) if j is odd. Consequently, we 
obtain from (3.2) (with p = 2d) 



where the last property follows from the orthogonality of the Legendre func- 
tions (see [1]). This implies (3.12) [or equivalently (3.10)] and proves the 
theorem. □ 

Remark 3.3. It should be noted that the mapping (8,<p) — > (sin#cos</>, 
sin 8 sin (p, cos 8) from the rectangle [0,ir] x (— tt,tt] onto the unit sphere S2 
maps all points of the form (0, (p) and (n, (p) with (p S (—71", 7r] onto the points 
(0,0,1) and (0,0,-1) on S%, respectively. Moreover, it is easy to see that 
the vector fd(8,(p) defined in (2.7) satisfies, for all cp € (— 7r,7r], 



As a consequence, various points of the designs \x ® v constructed by The- 
orem 3.2 can be identified on the unit sphere if the support of the factor \jl 
contains the point or n. To be precise, let p, denote the measure obtained 
from ji by omitting the points and tt, and define hq as the measure that 
puts masses ^({0}) and /j,({tt}) at the points (0,0) T and (tt,0) t , respectively. 
Then the measure 



has the same information matrix as the measure [i<®v. Note that in the case 
+ Ml 71 "}) = it follows that p, = fx, because the points and tt are 
not support points of the design \x. 

Example 3.4. Consider the spherical harmonic regression model of de- 
gree d = 1 and the case r = d + 1 = 2. From Lemma 3.1 with p = 2d = 2 it 
follows that the polynomial V%{x) = (x — l)(x + 1/3) satisfies 



J-i 

The points x\ = —1/3 and X2 = 1 generate a quadrature formular of degree 
2 on the interval [—1,1] with corresponding weights 




r 




/ d (oi) = (i,o,...,ofei( d+1 ) 2 , 
/ d (7r,0) = (-i,o,...,o) T e]R( d + 1 ) 2 . 



Po + A* ® v 
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^ arccos 

V; 





( 


7T 


7T 






~3 


3 


J O 




1 


1 


■) 




\ 


3 


3 


3/ 



According to Theorem 3.2 any design of the form 

, . / arccosf— i) \ . . 
//® i/(a,i) = I i 3 V 3J \(g)u(a,t) 

V 4 4 / 

with i > 3, a G (— ^7-71",— it] is ^-optimal for estimating the parameters 
K£ o k ^c (q < 1) and ^p.r-optimal in the first-order spherical harmonic re- 
gression model (2.6). A typical example is given by the six-point design 

I ^ 

4 4 

and by Remark 3.3 the design with equal masses at the points ( arccos (—|), — ^} 
(arccos(— |), ^), (0,0) and (arccos (— |), it) has the same information matrix 
as the design // (8) namely the identity matrix 14. 

Remark 3.5. Note that there are numerous possible ways to construct 
a discrete design with an information matrix equal to /(d+l) 2 m the spherical 
harmonic regression model (2.6). According to Theorem 3.2, a quadrature 
formula with r G {d + 1, . . . , 2(i} nodes si, . . . , x r with positive weights is 
required, which integrates polynomials up to degree 2d exactly. By (3.7), 
this formula gives the factor \i of the optimal design /i (8) u, where the second 
factor is any design of the form (3.8) with t > 2d + 1. By Lemma 3.1, the 
quadrature formula is determined by the equations 

(3.15) ^ V r (x)x e dx = 0, £ = 0,...,2d-r. 

If Pj(x) denotes the jth Legendre polynomial orthogonal with respect to 
the Lebesgue measure on the interval [—1,1], then the polynomial V r {x) 
can be represented as a linear combination of the Legendre polynomials 
Pq(x), . . . ,P r (x) and the orthogonality in (3.15) implies, for some constants 
a 2d-r+li • • • > a r> 

r 

(3.16) V r (x)= ajPjix). 

j=2d-r+l 

Note that the constants c^-r+i, ■ ■ ■ ,a r have to be chosen such that V r (x) 
has r real roots in the interval [—1,1] and such that the corresponding 
weights defined by (3.4) are positive. This is in general a nontrivial problem. 
However, one can easily describe a class of quadrature formulas for which this 
property is satisfied. For this, let P^ a ^'(x) denote the jth Jacobi polynomial 
orthogonal with respect to the measure (1 — x) a (l + xfdx on the interval 
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[—1,1] (see [22]). For any r > d + 1, it follows from these orthogonality 
properties that the identity (3.15) is satisfied for the polynomials 



Note that P r ' (x) is proportional to the Legendre polynomial P r {x) and 
consequently the representation of the form (3.16) is obvious, in this case 
choosing a r „i = • • • = a,2d-r+i = 0. Moreover, if p is the degree of one of these 
polynomials, it follows from classical results on orthogonal polynomials that 
each of the polynomials in (3.17) has precisely p roots in the interval [—1,1]. 
It is shown in [2] that the weights corresponding, by formula (3.4), to these 
roots are positive (see also [8]), and as a consequence we obtain a quadrature 
formula of degree 2d on the interval [—1,1], say {xj, 'Wj}j=i,...,p. The corre- 
sponding design fi obtained from (3.6) and (3.7) gives, in combination with 
any design v of the form (3.8), a <j? p -optimal design /i v for estimating the 
parameters k c in the spherical harmonic regression model (2.6). This 

design is also ^p^-optimal by Corollary 2.2. We finally note that the zeros 
of the polynomials in (3.17) are the support points of .D-optimal designs in 
heteroscedastic polynomial regression models (see [20]). 

Although there are numerous designs on the rectangle [0, ir] x (— tt, it] with 
information matrix in the spherical harmonic regression model (2.6) given by 
I(d+i) 2 -i the support points of the factor corresponding to the polar angle 9 
have to cover a sufficiently large range of the interval [0, it] in order to obtain 
a $ p -optimal design in the sense of Theorem 3.2. This is the statement of 
the final theorem of this section, which will be proven in the Appendix. 

Theorem 3.6. Let x\ denote the smallest zero of the Legendre poly- 
nomial Pd+i{x) and define z* = arccos If z > z* , then there exists no 
design on the rectangle [z,ir — z] x (— 7r,7r] with information matrix I[d+\) 2 
in the spherical harmonic regression model (2.6). 

4. Further discussion. 

4.1. The second-order spherical harmonic model. Consider the case d = 
2 corresponding to the second-order spherical harmonic model with nine 
regression functions. We calculate the designs corresponding to the four 
cases in (3.17) for the choice r = d + 1 = 3. From P^{x) = x(5x 2 — 3)/2 we 
obtain the support points of the probability measure \x corresponding to the 
polar angle 6 as 



(3.17) 



P^(x), (l-x)P^?(x), 

(i + x)pf^ (x), (i - x 2 )p r y } (x). 




02 = arccos = — 
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Fig. 3. Projections of the suppport points of the optimal design in the second-order spher- 
ical harmonic regression model (2.6). The left panel corresponds to the optimal design 
defined by (4.1), while the right panel represents the design (4.3). 



and the corresponding weights are given by 

f A o\ 1 f 1 ^ 2 -3/5 , 4 l-w 2 5 

(4.2) W2 = - —, — ax = —, wi=W3~ 



2j-i (-3/5) 9 2 18 

Similarly, if the polynomial (x 2 — tyPy^ (x) = j(x 2 — l)(5x 2 — 1) is used for 
the construction of a quadrature formula, we obtain 



#1 = arccos 1 = 0, 62= arccos \ — , 

(4-3) V5 

#3 = arccos ( — J , 64 = arccos(— 1) = tt, 

and the weights are obtained from the representation (3.4) and given by 
(4.4) w 1 = W 4 = -L, W2 = W3 = JL. 

The measure v corresponding to the azimuthal angle (j> is given by (3.8), 
where t > 5. The projections of the support points of the two measures onto 
the unit sphere are depicted by the left and right panels of Figure 3, where 
for the second component a design with t = 5 support points and a = —tt is 
used. It should be noted that it follows from Remark 3.3 that the design //(g) v 
obtained from (4.3) and (4.4) for the factor [i and v = v(a, 5) corresponds to 
uniform design on the sphere with 12 support points. A general construction 
of uniform designs will be discussed in the following paragraph. We finally 
illustrate the two other non-symmetric cases in (3.17). For the supporting 
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polynomial we obtain 

(1 - x)P 2 (1 ' 0) (x) = ^y^($x 2 + 2x - 1), 
which gives, for the support points of the factor [i of the product design 

6 1 = arccos 1 = 0, 

l + v^ 



(4.5) 02 = arccos 



arccos 



5 



with corresponding weights 

101 = 5, 102 = 3^(16 + ^6), 103 = ^(16-^). 

The fourth case of a factor ^ of a design v(a, t) yielding Ig as information 
matrix in the second-order spherical harmonic model (2.6) is obtained by 
symmetry, that is, 



arccos 



l + vg 

5 



(4.6) 6>2 = arccos^ — ^r—^J, @3 = arccos(l) = it, 

wi = ^(16- V6), ^2 = ^(16 + ^), ^3 = 3- 

The projections of the support points onto the unit sphere of the corre- 
sponding product designs /x x v(a, t) with a = —ir and t = 5 are depicted by 
the left and right panels of Figure 4. Note that the two cases are related by 
a reflection of the support points at the equator. 



4.2. Optimal designs with equal weights. Note that the designs provided 
by Theorem 3.2 are in general not uniform designs, which would have equal 
weights at their support points. Because designs with this structure are par- 
ticularly attractive from a practical point of view, we will briefly discuss the 
possibility of their construction in this section. Note that for the determina- 
tion of an optimal design in the spherical harmonic regression model (2.6) 
with equal weights at its supports by the procedure introduced in Section 3, 
it is necessary to find a quadrature formula that has equal weights at its sup- 
port points and integrates polynomials of degree 2d exactly. This problem 
has a long history in mathematics (see [6]). 

It follows from Lemma 3.1 that such a formula must have at least d + 1 
nodes. Moreover, quadrature formulas with the minimal number of n nodes 
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Fig. 4. Projections of the suppport points of the optimal design in the second-order spher- 
ical harmonic regression model (2.6). The left panel corresponds to the optimal design 
defined by (4.5), while the right panel represents the design (4.6). 

and equal weights which integrate polynomials of degree n exactly only exist 
in the cases n = 0, 1, . . . , 7, 9 (see [6], page 58). In the present context these 
formulas correspond to uniform designs for the polar angle 9 in the case 
d= 1,2,3,4, which give, in combination with a design of the form (3.8), a 
or ^p r-optimal design in the spherical harmonic regression model (2.6) 
with equal weights at its support points. The corresponding nodes of the 
quadrature formula required for the factor [i are depicted in the first four 
rows of Table 1. If d > 5, quadrature formulas with equal weights on 2d 
points integrating polynomials of degree 2d exactly do not exist and more 
nodes are required for the construction of such formulas. We determined 
formulas of this type numerically for d = 5, 6, 7 and depicted them in the 
last three rows of Table 1. Note that these formulas use the origin and 
that the number of nodes increases rapidly. For example, if d = 7, we only 
found a quadrature formula with 23 support points and equal weights which 
integrates polynomials of degree 14 exactly. 

4.3. Some efficiency considerations. It might be of interest to compare 
some of the commonly used designs from the literature on the analysis of 
3D shapes with the ^ p , r - an d ^-optimal designs obtained in the present 
paper. As an example, we consider two uniform designs of the form jjl® u, 
where the components are given by 



(4.7) 
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Table 1 

Quadrature formulas with equal weights at their nodes* 

d X ±CCi ±X 2 ±X 3 ±X4 ±IE5 ±X G ±X7 ±X S ixg ±CCio ±X 



11 



1 0.577 

2 0.188 0.795 

3 0.267 0.423 0.866 

4 0.168 0.529 0.601 0.912 

5 0.223 0.247 0.443 0.671 0.724 0.939 

6 0.008 0.282 0.358 0.458 0.566 0.760 0.778 0.954 

7 0.174 0.177 0.186 0.328 0.502 0.533 0.542 0.712 0.797 0.852 0.965 

*These formulas correspond to a uniform design jj, u of the form (3.7) for the polar angle 
6. The resulting design /n u (g) v(a,t) with v(a,t) defined by (3.8), t > 2d + l, is a 4> p - or 
^p.r-optimal design in the spherical harmonic regression model (2.6) of degree d with 
equal weights at its support points. 

and the designs are given by 

2i \\ / 2z7r 



arccos 1 



(4.8) fj, ■ 




/ ZITT 
I 7T 

1 



i=l,...,7ll 1^2 / i=l,...,Tl2 



respectively. Note that the design // (g> v defined by (4.7) corresponds to a 
uniform distribution on the grid in the rectangle [0, 7r] x [— 7r,7r], while a 
design \t% v of the form (4.8) yields a design on the sphere that takes obser- 
vations on several circles with equal distance on the z axis. This design was 
used by Ding, Nesumi, Takano and Ukai [5] for a principal component anal- 



Table 2 

D, E, A and 9-i >r efficiencies (r = 2,3) of the uniform designs (4.7) and (4.8) in first- 
and second-order spherical harmonic models 



Design (4.7) Design (4.8) 



d 


m 


eff D 


eff B 


e« A 


eff*_ 12 


eff *-i, 3 


eff D 




e« A 


eff*_ 12 


eff *-l,3 


1 


3 


1.000 


1.000 


1.000 


1.000 


1.000 


0.940 


0.500 


0.870 


0.667 


0.789 




4 


0.997 


0.938 


0.994 


0.938 


0.957 


0.964 


0.600 


0.923 


0.750 


0.857 




5 


0.993 


0.900 


0.986 


0.900 


0.931 


0.976 


0.667 


0.949 


0.800 


0.894 




6 


0.989 


0.875 


0.979 


0.875 


0.913 


0.983 


0.714 


0.964 


0.833 


0.916 




7 


0.986 


0.857 


0.973 


0.857 


0.900 


0.987 


0.750 


0.973 


0.857 


0.931 


2 


4 


0.991 


0.801 


0.982 


0.838 


0.851 


0.902 


0.229 


0.745 


0.331 


0.427 




5 


0.987 


0.805 


0.974 


0.824 


0.831 


0.935 


0.323 


0.838 


0.435 


0.539 




6 


0.981 


0.794 


0.964 


0.807 


0.811 


0.954 


0.399 


0.888 


0.512 


0.617 




7 


0.976 


0.782 


0.955 


0.793 


0.796 


0.965 


0.461 


0.918 


0.571 


0.674 




8 


0.972 


0.772 


0.947 


0.782 


0.785 


0.973 


0.511 


0.937 


0.617 


0.716 
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Table 3 

D, E, A and *I/_i,r efficiencies (r = 2,3) of the designs (4.7) and (4.8) in spherical 
harmonic models of degree 3 and 4 



Design (4.7) Design (4.8) 



d 


ni 


eff D 


eff B 


e« A 


eff*_ lj2 


eff*_ 1 3 


eff D 


effB 


effA 


eff*_ 12 


eff*_ 1 3 


3 


5 


0.980 


0.799 


0.961 


0.802 


0.803 


0.874 


0.094 


0.600 


0.146 


0.199 




6 


0.975 


0.784 


0.953 


0.784 


0.787 


0.911 


0.155 


0.733 


0.223 


0.296 




7 


0.970 


0.768 


0.944 


0.768 


0.772 


0.934 


0.214 


0.811 


0.292 


0.377 




8 


0.965 


0.756 


0.936 


0.756 


0.760 


0.948 


0.269 


0.859 


0.352 


0.444 




9 


0.961 


0.747 


0.929 


0.747 


0.751 


0.959 


0.318 


0.891 


0.404 


0.500 


4 


6 


0.969 


0.739 


0.942 


0.755 


0.761 


0.851 


0.035 


0.434 


0.057 


0.081 




7 


0.965 


0.747 


0.936 


0.751 


0.753 


0.890 


0.067 


0.600 


0.102 


0.142 




8 


0.961 


0.739 


0.929 


0.742 


0.743 


0.916 


0.103 


0.709 


0.149 


0.203 




9 


0.956 


0.731 


0.922 


0.733 


0.734 


0.933 


0.142 


0.781 


0.196 


0.261 




10 


0.952 


0.724 


0.915 


0.726 


0.727 


0.945 


0.180 


0.830 


0.240 


0.314 



ysis of the variance-covariance matrix of the least squares estimator (1.3) in 
the spherical harmonic regression model of degree 7. In Tables 2 and 3 we 
consider the problem of designing an experiment to estimate the full param- 
eter vector c or the principal component analysis in the spherical harmonic 
regression model (2.6). We show the D, E and A efficiencies of these designs 
for various values of n± and ri2 in the spherical harmonic models of degree 1 , 
2, 3 and 4. The tables also contain the ^_i jr -efficiencies, which are defined 
by 

M) Bff ff|P , r (fl= % f\ y 

For the sake of brevity we choose the design for the azimuthal angle (j> 
as the design defined by (3.8) — other uniform designs for this component 
will yield substantially lower efficiencies and are therefore not depicted. As 
a consequence, the efficiencies of the design /i (g> v do not depend on rti 
(provided that ni > 2d-\- 1, which will be assumed throughout this section). 

For the first-order spherical harmonic regression model, we observe very 
good D and A efficiencies of the designs (4.7) and (4.8). However, the E 
efficiencies and the ^ Pt r efficiencies of these designs are substantially smaller, 
in particular those efficiencies obtained for the design (4.8) with moderate 
values of n\. For spherical harmonic models of larger degree both designs will 
still yield high D efficiencies, the designs defined by (4.7) yield reasonable 
A efficiencies, but the E and V^—i^ efficiencies are substantially smaller. 
Moreover, the A and E efficiencies of the design (4.8) are very low. The Vf-i.r 
efficiencies of this design are slightly larger but still not satisfactory. It is 
also interesting to note that the efficiencies of the design (4.7) are decreasing 
with the number n\ while they are increasing for the design (4.8). 
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From these calculations and additional results, which are not depicted for 
the sake of brevity, we observe that designs with a uniform distribution on 
a grid in the rectangle [0,tt] x (— 7r,7r] should be used only if maximization 
of the D criterion is the preliminary goal of the design of the experiment. 
Whenever principal component analysis is the main goal of the experiment 
or precise estimates of the parameters themselves are required, some more 
care is necessary in the design of an experiment for analyzing 3D shapes. In 
this case, the loss of efficiency when using uniform designs on a grid in the 
rectangle [0, tt] x [— ir, tt] or a uniform design taking observations on several 
circles with equal distance on the z axis can be substantial. In most cases 
there exist substantially more efficient designs for the analysis in a spherical 
harmonic regression model. The advantages of the optimal designs derived 
in the present paper will also be illustrated in the following example. 

4.4. A concluding example. To demonstrate the benefits of our designs, 
we finally reanalyze the design used by Ding, Nesumi, Takano and Ukai [5] 
for shape analysis of Citrus species. These authors used observations at 360 
points measured in 10 circles using the equal height sampling method for 
the z axis. By choosing d = 7 the data for the surface shape were expanded 
into the first 64 terms of spherical harmonic functions. The information con- 
tained in the spherical harmonic coefficients was summarized by a principal 
component analysis using the first seven principal components. Note that the 
design of Ding, Nesumi, Takano and Ukai [5] corresponds to a design of the 
form (4.8) with n\ = 10 and n2 = 36. In the first row of Table 4 we show the 
efficiencies of this design with respect to the optimal designs obtained in this 
paper. The first factor of an optimal (uniform) design can be obtained from 
the quadrature formula corresponding to the spherical harmonic regression 
model of degree 7 in Table 1. We observe a reasonable efficiency only for the 
.D-criterion. For all other criteria the design used by Ding, Nesumi, Takano 
and Ukai [5] is very inefficient. The optimal designs proposed in this paper 
(or appropriate approximations) will yield substantially smaller variances of 
the least squares estimates for linear combinations of the parameters. 

We finally note that the optimal designs obtained in this paper are ap- 
proximate and do not have masses that are multiples of 1/360. However, 
a very efficient design for the inference in the spherical harmonic regres- 
sion model of degree 7 using 360 different points can easily be obtained as 
follows. The quadrature formula obtained from Table 1 has equal masses 
at 23 nodes and yields a design of the form (3.7), say {^,1/23}?^, where 
9i = arccosxj and xi,...,X23 are the points listed in Table 1 for the case 
d = 7. We propose to combine this design with two designs of the form (3.8) 
to obtain an efficient exact design. More precisely, we propose to use the 
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Table 4 

D, A, E and efficiencies eff*^ r (r = 1, . . . , 10) in the spherical harmonic regression 
model of degree 7 used by Ding, Nesumi, Takano and Ukai [5]* 













eff ^ 


-l,r 












1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


eff A 


eff D 


0.003 
0.958 


0.006 
0.958 


0.008 
0.958 


0.011 
0.958 


0.013 
0.958 


0.016 
0.958 


0.019 
0.958 


0.021 
0.958 


0.024 
0.958 


0.026 
0.958 


0.149 
0.987 


0.840 
0.992 



*The first row is the design given by (4.8) with n 1 = 10 and = 36. The second row gives 
the efficiencies of the design defined by (4.10), which takes one observation at each point 
of the set IA. 



uniform distribution at the points 
(4.10) 



U 



15 

vr(2j - 16) 



16 



i = l,...,8;j = 1, . . . , 15 
i = 9,...,23;j = l,...,16 



as the design for the spherical harmonic regression model of degree 7. The 
implementation of this design is illustrated in Figure 5. The efficiencies of 
the exact design £jy are depicted in the second row of Table 4 and we observe 
that this design, which advises the experimenter to take one observation at 
each point of IA, is highly efficient with respect to all optimality criteria 
under consideration. 



5. Conclusions and some directions for future research. In this paper 
we tried to provide a starting point for studying design problems for 3D 
shape analysis. There are many issues in this area which can be addressed 
by the choice of an experimental design and we concentrated on spherical 
harmonic descriptors, with special emphasis on classification. This direction 
was motivated by recent work of Ding, Nesumi, Takano and Ukai [5], who 
used the coefficients in a spherical harmonic expansion to classify different 
fruit shapes by principal component analysis. Other authors compare dif- 
ferent shapes by simply computing the Euclidean distance between their 
corresponding spherical harmonic descriptors (see [7]). If classification is the 
main object in 3D shape analysis, precise estimation of the coefficients in 
the spherical harmonic expansion is of particular importance and should be 
addressed by the choice of an experimental design, whenever this is possible. 
In the present paper we have constructed optimal designs for this purpose by 
considering Kiefer's $ p criteria and a new criterion which is directly related 
to principal component analysis. It is demonstrated that the new designs 
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Fig. 5. Implementation of the design (4.10) on the unit sphere. 

are substantially more efficient for parameter estimation and principal com- 
ponent analysis than the commonly used designs by reanalyzing a typical 
data example from the literature. 

As pointed out by a referee, there are several other design issues concern- 
ing image analysis that are not addressed in this work, but are important 
topics for future research in this direction. In the following list we briefly 
mention the — in our opinion — most important directions in image analysis, 
where the generalization of our approach could be profitable. 

1. The main motivation of our work stems from the fact that spherical 
harmonic descriptors are used for classification by principal component 
analysis. There are many other statistical tools for this purpose, for ex- 
ample, independent component analysis (see the recent monograph by 
Hyvarinen, Karhunen and Oja [10]). Therefore, it is important to deter- 
mine efficient designs for alternative classification rules and to compare 
these with the designs derived in this paper. 

2. In addition to the problem of classification, another important goal of 
3D shape analysis is the reconstruction of shapes. In such cases the num- 
ber of basis functions is usually very large and an optimality criterion 
should focus on the closeness between the original shape and the re- 
constructed shape. A suitable optimality criterion for this purpose is to 
minimize the expected mean squared error with respect to the choice of 
an experimental design. This problem is closely related to the question 
of model uncertainty. Since the pioneering work of Box and Draper [3], 
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several authors have addressed the problem of incorporating the bias in 
the optimality criterion (see, e.g., Wiens [24] and the references therein). 
We note, however, that by Theorem 2.1 the uniform distribution on the 
sphere is <£ p -optimal in the model (2.6) for any order d 6 N. Therefore, 
from a theoretical viewpoint, the question of model uncertainty does not 
exist for the uniform distribution on the sphere (it is the optimal design 
for any d £ N). However, these problems arise if a design has to be im- 
plemented for a fixed sample size. If the main goal of the analysis is the 
reconstruction of shapes, we recommend choosing d as large as possible. 
In Section 4.4 we have indicated how such a design can be found in the 
case d = 7, which was given in a concrete application of spherical har- 
monic descriptors in biology. In principle this method can be adapted to 
larger values of d. However, the numerical difficulties of finding quadra- 
ture formulas of degree 2d with equal weights increase substantially with 
the value d and suitable software has to be developed for this purpose. 
We also note that the sample size n always provides an upper bound for 
the degree of the model (2.6). Therefore — in principle — the problem of 
model uncertainty is caused by a limited sample size, not by our approach 
of constructing optimal designs. 

3. It is remarkable that the approach developed in this paper does not re- 
quire any distributional assumptions, because it is based on the least 
squares method and the covariance matrix of the corresponding estima- 
tor. However, the assumption of uncorrelated observations is important 
for the calculation of this matrix. Observations from image analysis are 
usually taken from one object and therefore this assumption has either to 
be checked or to be taken into account by the construction of the design. 
Note that many authors working in image analysis do not consider an er- 
ror in distribution. In our model (2.6), the error refers to a nonsystemtatic 
measurement error, for which the assumption of no correlation may be 
justified. In general, the construction of optimal designs for models with 
correlated observations is a very hard problem even in less sophisticated 
models than considered in this paper (see, e.g., [18]). The incorporation 
of correlation structures in the construction of optimal designs is proba- 
bly one of the most challenging problems in future research on optimal 
designs for image analysis. 

4. It seems to be desirable to extend the methods derived in this paper to 
other problems in image analysis, in particular to the problem of 2D image 
reconstruction. We note that the design problem for 3D shape analysis 
is substantially simpler than design problems for general image analysis. 
The reason for this is the specific structure of the system of spherical 
harmonic functions used in (1.1), which is not available for 2D image 
analysis. A reasonable system for the two-dimensional case is the Zernike 
polynomials (see [14]), and research to construct good designs in this 



26 



H. DETTE, V. B. MELAS AND A. PEPELYSHEV 



case is still ongoing. Initial results in this direction indicate that there 
are substantial differences between the two- and the three-dimensional 



cases. 



APPENDIX: TECHNICAL DETAILS 

Proof of Lemma 3.1. Assume that conditions (A) and (B) of Lemma 3.1 
are satisfied and let Q(x) denote an arbitrary polynomial of degree z. By 
Bezout's theorem the polynomial Q can be represented in the form 

(A.l) Q(x) = P(x)V r (x) + R(x), 

where V r (x) = JYj=i( x ~ x j)i the polynomial P{x) is of degree z — r and the 
degree of R(x) is less than r. Because the degree of R is at most r — 1, it 
can be represented as 

r 

(A.2) R(x) = J2^(x)R(xj), 

3=1 

and we obtain from the conditions (A) and (B) of the lemma that 
\ J Q(x)dx = ^ J R(x)dx 



.7=1 J ~ 1 



= S R ( x j) w J = Yl Q( x i) w r 

3=1 3=1 

Using the functions Q(x) = x e (£ = 0, . . . ,z — r) yields the identities in (3.2). 

For a proof of the converse, we assume that (3.2) is valid and obtain for 
£ = 0,...,z-r, 



f 1 r 
/ V r (x)x dx = ^ V r (xj 



)xjWj = 0, 

3=1 

which gives condition (A). On the other hand, condition (B) follows from 
the identity 



(Xi) = Wj, 



observing the property ij(xi) = 5{j for the Lagrange interpolation polynomi- 
als. □ 
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Proof of Theorem 3.6. Let £ denote an arbitrary discrete design 
on the rectangle [0, 7r] x (—it, it] and denote by 9\,...,9n the distinct first 
coordinates of the corresponding support points of £. Obviously £ can be 
represented as 

N 

(A.3) £ = 

i=l 



6 



where the designs & are defined by 

h,<t>il) ••• (9i,<friNi 

wn ■■■ w iNi 

(with iVj E N) and represent the part of £ corresponding to the support 
points with first coordinate equal to 9i (i = 1, . . . , iV). Define Xi = cos^j and 
lUi = J2j=i w ij (« = !,•••, -/V), and consider the design 



-.7 

(x\ ■■■ X N 
\Wl ■■■ W N 



(A.4) V = Tfe 



If the design £ satisfies the condition M(£) = -f(d+i) 2 j it follows for the 
submatrix corresponding to the (d + 1) regression functions Y®(9,(f)) = 
^£ + Tif (cos 0) (£ = 0, . . . , d) that 



d 

I 1 / Y?(0,<l>)Y£(6,<l>)d£(0,<l>) 

V0 ./-• 7r / i,k=0 

/ Af AT* \ d 

(A.5) = ^^^V27TTv / 2fcTlP^ (x i )P fc (x l ) 

\ i=l j=l / £,/e=0 

/ N \ d 

= V2l+TV2k + TJ2mPi(xi)P k (xi) 

V i=l / £,k=0 

[Note that Pf(x) = Pi(x) is the ith. Legendre polynomial.] On the other 
hand, the orthogonality relation for the Legendre polynomials (see [22]) 
yields 

(A.6) }-j\( x ) Pk ( x ) dx = J^ i £,k = 0,...,d, 

and it follows from (A.5) that the quadrature formula defined by the de- 
sign n in (A.4) integrates polynomials of degree Id exactly. The assertion of 
Theorem 3.6 now follows from the following auxiliary lemma. □ 



Lemma A.l. For any n £ N, let E denote the set of all probability mea- 
sures on the interval [—1,1] of the form (AA) which integrate polynomials 
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of degree n exactly, that is, 



(A.7) i / x l dx = Y J WjX% i = 0, 



Then 



minmin{t> G R + | suppr/ C [— v,v]} = \x*\, 



where x\ is the smallest zero of the rth Legendre polynomial P r {x) and r 
|_n/2j +1. Moreover, the minimum is attained by the measure 



(A.8) n* 



X ^ X y, 

Wi •■• W„ 



where the points the zeros of the Legendre polynomial P r {x), 

Wi = 7j /_ 1 ^j(x) (i = 1, . . . , r) and is £/ie Lagrange interpolation poly- 
nomial with nodes x*,...,x*. 

Proof. For a design 77 on the interval [—1,1], let 

(A. 9) 14(77) := nrin{t> £ M + | supp?? C [— v, t>]} 

denote the half of the minimal length of intervals that contain the support 
of 77 and define for each JVgN the set E at C 3 as the set of all probability 
measures with N support points satisfying (A.7) (note that E = Ujvg^Sjv). 
Because for a design 77 of the form (A. 4) the equations in (A.7) can be 
written as 



(A.IO) 



it follows from Caratheodory's theorem (see [19]) that there exists a measure, 
say 77, with at most n + 1 support points Xj 1 , . . . , Xj n+1 (1 < j\ < ji < ■ ■ ■ < 
jn+i < N) such that (A.7) is satisfied. Consequently, we obtain 

(A. 11) inf 7^(77) = inf u{rj). 

If 

„ - f Xl " ' Xn + 1 \ a n 

V^i ••• Wn+iy 
with — 1 < x\ < • ■ • < x n+ \ < 1 is any probability measure on the inter- 
val [—1,1] that satisfies (A.7), then by Lemma 3.1 the weights can be rep- 
resented as 




Wi = \ i £i(x)dx, 
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where li(x) is the ith Lagrange interpolation polynomial with nodes x±, . . . , 
x n -\-\. Assume that |xi| = u{rj) and, for e > 0, define fj as the measure with 
weights 

^ .± n+l ^ 

(A.I2) wt = - ^—^dx, i = l,...,n + l, 

L J — 1 . ... Xi Xj 

at the points x\ = x± + e, xi = X2, ■ ■ ■ , x n+ \ = x n+ \ . If e is sufficiently small, it 
follows that all weights Wj are positive, which implies fj £ H n+ i and xi > x\. 
In the case |x n +i| = u{rf) we apply exactly the same argument to the point 
x n+ \ and obtain a measure fj with 14(77) < n ( 7 ?)- Consequently, the infimum 
in (A. 11) cannot be attained in H n+ i, that is, 

(A. 13) inf 7/(77) = i n f -u ( 7 ?)- 

We now prove that the infimum on the right-hand side cannot be attained 
in 3t whenever t > [n/2\ + 1 = r. For this, consider a measure 

(x\ ■■■ x t 
\W\ ••• w t 



Recall from Lemma 3.1 that the weights have to satisfy (3.4) (with t = r) 
and that by (3.3) the support points satisfy 



»i t 

(A.14) / l(x-x j )x i dx = 0, 



0, . . . ,n — t. 



Without loss of generality we assume that \x±\ = 74(77) and we will again 
construct a design with a smaller value of 74(77). For this we define x\ =x±+e, 
Xj = Xj ( j = 3 + n — t , . . . , t) , 

t 

(A. 15) tp(x) = (x — xi) (x — ij) 

j=3+n—t 

and construct the remaining points, say X2, ■ ■ ■ , x n -t+2, such that the result- 
ing measure with weights of the form (3.4) at the points x\, . . . ,xt defines a 
quadrature formula of degree n. A necessary condition for this property is 
given by 

■1 n-t+2 

(A.16) / ip(x) TT (x- Xj )x e dx = 0, £ = 0,...,n-t, 

J-l j=2 

and a straightforward calculation shows that these equations are equivalent 
to the system 

Q r\ n-t+2 

(A.17) — / if>(x) J| (x-xj) 2 dx = 0, i = 2,..., n-t + 2, 

j=2 
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which determines the points Xj = Xj(e) (j = 2, . . . , n — t + 2) implicitly as a 
function of the parameter e, where Xj(0) = Xj (j = 2, . . . ,n — t + 2). We now 
consider the Jacobi matrix of the system (A. 17), 

/ Q2 ,1 n-t+2 \ n-t+2 

(A.18) J(x,e)=(—— ^x)Y[(x-x k ) 2 dx) 

yOXjOXjJ-! k=2 J.. =2 

where we use the notation x = (£2, ■ ■ ■ ,x n -t+2)- This matrix can be calcu- 
lated using the quadrature formula corresponding to the measure r\. For this, 
note that x(0) = (x2, ■ ■ ■ ,x n _t+2) and define the polynomial 

n-t+2 

(A.19) gi(x) = 2ip(x) Y[ {x-Xjfdx. 

Note that the degree of gi is less than or equal to n. If a € W denotes the 
ith. unit vector, we obtain 

ef J(x(0),0)ei = / g i (x)dx = y2g i (xj)wj = g i (x i )w i 

t n-t+2 
= 2w i (x i -X 1 ) Y[ (Xi-Xj) Y[ (Xj-Xif^O, 
j =n -t+3 j=2,jj^i 

and by a similar calculation, e[ J(x(0), 0)ej = whenever i 7^ j (i, j = 2, . . . , 
n — t + 2). Consequently, J(x(0),0) is diagonal with det J(x(0), 0) / 0. It 
now follows from the implicit function theorem (see [9]) that for sufficiently 
small e there exist analytic functions x^^s), . . . ,x n —t+2{£) such that (A. 16) 
is satisfied for the points x\,...,xt- This implies that for sufficiently small 
e, the weights defined by (A. 12) are positive and the design 

( x\ ■■■ x t 
Ui • • • Wt. 



(A.20) fj 



defines a quadrature formula of degree n with positive weights, that is, 
fj E 3$. If necessary [in the case \xt\ = u{rf)\ we use a similar argument for 
the largest support point of the probability measure r] and finally obtain 
a probability measure fj E St such that u(fj) < u(jj). In other words, the 
infimum on the right-hand side of (A. 13) cannot be obtained in if t > 
[n/2j +1. From (A. 15) it is easy to see that this construction is not possible if 
t < r = [ra/2j + 1. Moreover, it is well known (see, e.g., [8]) that a quadrature 
formula that integrates polynomials of degree n exactly must have at least r 
support points. Moreover, if it has r support points, it is uniquely determined 
and given by the probability measure rf defined in (A. 8), which completes 
the proof of Lemma A.l. □ 
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